19  Production model - river 4 (IS)

Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.

In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv

Model code is in ./models/production/modelProduction_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/modelProduction_OB_makeFile.R
Functions for this qmd file in ./models/qmdProduction_OB_functionsToSource.R

There are 3 models here: All cohorts:
model 1, with quadratic weight effect on phi By cohort:
model 2, with linear weight effect on phi model 3, with quadratic weight effect on phi

19.0.1 How many ageInSamples to include?

all <- tar_read(cdWB_electro_target)
#str(all)
table(all %>% filter(river == "wb obear")|> dplyr::select(ageInSamples))

ageInSamples
   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
1308 1890 1007  837 1486  875  446  298  444  159  155   57   47   17   16    3 
  16   18 
   3    3 

table(all %>% filter(river == "wb obear")|> dplyr::select(cohort,ageInSamples))

      ageInSamples
cohort   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  18
  1999   0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   0
  2000   0   0   0   0   0   0   0   0  19  10   0   6   3   1   1   1   1   0
  2001   0   0   0   0  79  73   0  58  48  22  14  15   8   3   3   1   1   3
  2002  73 131   0 120 137  83  35  59  41  13  33   6   3   1   7   1   1   0
  2003 104 147  52 118 129  86  49  50  56  22  73   5   6   4   3   0   0   0
  2004   3  95   0  44 109  62 100  18  28  17   8   4   6   4   1   0   0   0
  2005  46  54  42  18  31  32  14   5   4   1   0   0   0   0   1   0   0   0
  2006  21  73  61  27  60  33  16   9   9   5   2   0   0   0   0   0   0   0
  2007  28  49  46  26  52  28  25   5   6   3   2   2   0   0   0   0   0   0
  2008  11  44  52  21  28  14  18   8  14   1   2   1   1   0   0   0   0   0
  2009  44 310 213 154 225 122  84  38  63  18  11  14   5   0   0   0   0   0
  2010  89  36  28  13  21  16  12  11   5   2   0   0   0   0   0   0   0   0
  2011   8  21  22  24   9  11   8   4   5   0   0   0   0   0   0   0   0   0
  2012 167 307 319 197 196 142  72  28  39  20  10   4   4   3   0   0   0   0
  2013  17  67  63  18  33  24  13   5  10   3   0   0   0   0   0   0   0   0
  2014  31 168 109  57 100  73   0   0  17   0   0   0   0   0   0   0   0   0
  2015 395 295   0   0 113   0   0   0   1   0   0   0   0   0   0   0   0   0
  2016 157   0   0   0  56   0   0   0   4   5   0   0   6   0   0   0   0   0
  2017   6   0   0   0  28  36   0   0  10   0   0   0   3   0   0   0   0   0
  2018  45  16   0   0  21   0   0   0  53   2   0   0   1   1   0   0   0   0
  2019  51   0   0   0   2  19   0   0  12  15   0   0   0   0   0   0   0   0
  2020  12  44   0   0  57  21   0   0   0   0   0   0   0   0   0   0   0   0
  2021   0  33   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0

cohorts <- 2002:2014

19.1 All Cohorts, modelNum 1: phi(i,t,c) * g(i,t,c), p(i,t) model

19.1.1 Model code

# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_1_mostRecent.RData') # WILL NEED TO CHANGE TO MOST RECENT
#load('./models/production/runsOut/production_OB_model_2_2025-02-28 19.RData')

d$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] <- grIntT[t] + grBeta[1, cohort[i]]
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 0.5)
    }
    for (c in 1:nCohorts) {
        grBeta[1, c] ~ dnorm(0, sd = 1/0.667)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t] + phiBeta[2, i, t] * weight[i, 
                t] * weight[i, t] + phiBetaCohort[cohort[i]]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:2) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (c in 1:nCohorts) {
        phiBetaCohort[c] ~ dnorm(0, sd = 1/0.667)
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
        phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        for (t in 1:(T - 1)) {
            weightZ01[i, t] <- weightIfAvailable(weight[i, t], 
                sdWeight, meanWeight, z[i, t], available[i, t])
            z01[i, t] <- aliveIfAvailable(z[i, t], available[i, 
                t])
        }
    }
    for (t in 1:(T - 1)) {
        sumWeightZ01[t] <- sum(weightZ01[1:N, t])
        sumZ01[t] <- sum(z01[1:N, t])
    }
    for (cohortLoop in 1:nCohorts) {
        for (t in 1:(T - 1)) {
            sumWeightZ01ByCohort[cohortLoop, t] <- sumVecForCohort(weightZ01[1:N, 
                t], cohortLoop, cohort[1:N])
            sumZ01ByCohort[cohortLoop, t] <- sumVecForCohort(z01[1:N, 
                t], cohortLoop, cohort[1:N])
        }
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

d$runTime

Time difference of 6.697066 days

  MCMCplot(object = d$mcmc, params = c("sumZ01"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("sumWeightZ01"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("sumZ01ByCohort"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("sumWeightZ01ByCohort"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("p"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("grIntT"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("phiIntTOut"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("phiBetaT"), ref_ovl = TRUE)

  MCMCplot(object = d$mcmc, params = c("phiBetaCohort"), ref_ovl = TRUE)

MCMCsummary(object = d$mcmc, params = c(
        #"phiIntT", 
        "phiIntTOut", 
        "p", 
        "grIntT", 
        "phiBetaT",
        "sumWeightZ01",
        "sumZ01"
      ), round = 3) %>%
      rownames_to_column(var = "var")

                var      mean      sd      2.5%       50%     97.5%  Rhat n.eff
1     phiIntTOut[1]     0.751   0.075     0.600     0.757     0.899  1.13    33
2     phiIntTOut[2]     0.707   0.090     0.513     0.712     0.865  1.22    19
3     phiIntTOut[3]     0.502   0.104     0.321     0.502     0.687  3.57    31
4     phiIntTOut[4]     0.378   0.073     0.242     0.374     0.515  1.99    25
5     phiIntTOut[5]     0.522   0.091     0.334     0.536     0.684  2.19    22
6     phiIntTOut[6]     0.579   0.096     0.381     0.577     0.765  1.89    18
7     phiIntTOut[7]     0.743   0.061     0.607     0.743     0.850  1.01    22
8     phiIntTOut[8]     0.465   0.105     0.300     0.449     0.686  1.00    12
9     phiIntTOut[9]     0.452   0.135     0.156     0.470     0.681  3.42    16
10   phiIntTOut[10]     0.469   0.115     0.269     0.477     0.671  2.75    16
11   phiIntTOut[11]     0.527   0.128     0.277     0.513     0.805  3.13    14
12             p[1]     0.597   0.023     0.553     0.597     0.641  1.22   217
13             p[2]     0.433   0.019     0.396     0.432     0.471  1.25   183
14             p[3]     0.708   0.019     0.670     0.708     0.744  1.01   338
15             p[4]     0.655   0.023     0.611     0.655     0.697  1.00   183
16             p[5]     0.627   0.026     0.575     0.627     0.679  1.03   130
17             p[6]     0.500   0.028     0.446     0.500     0.556  1.12   290
18             p[7]     0.709   0.037     0.636     0.711     0.779  1.20    79
19             p[8]     0.787   0.047     0.686     0.789     0.873  1.01   311
20             p[9]     0.605   0.068     0.468     0.606     0.731  1.01   182
21            p[10]     0.692   0.082     0.526     0.695     0.839  1.19   210
22            p[11]     0.835   0.112     0.591     0.850     0.994  1.04    38
23        grIntT[1]     0.158   0.088     0.045     0.155     0.275 13.83    76
24        grIntT[2]     0.201   0.084     0.096     0.196     0.309 16.30    59
25        grIntT[3]     0.936   0.085     0.824     0.930     1.051 12.12    73
26        grIntT[4]     0.392   0.084     0.285     0.382     0.498 15.50    60
27        grIntT[5]     0.153   0.086     0.038     0.154     0.271 11.14   155
28        grIntT[6]     0.277   0.084     0.168     0.285     0.389 11.65    47
29        grIntT[7]     1.677   0.087     1.554     1.667     1.803  9.05   187
30        grIntT[8]     0.791   0.088     0.668     0.795     0.919  9.34   133
31        grIntT[9]    -0.147   0.091    -0.300    -0.146     0.002  4.92   396
32       grIntT[10]     0.472   0.103     0.292     0.476     0.640  4.21    91
33       grIntT[11]     2.262   0.155     1.945     2.268     2.536  1.43    85
34   phiBetaT[1, 1]     0.044   0.602    -1.054     0.008     1.213  1.02    13
35   phiBetaT[2, 1]     0.929   0.671    -0.735     1.011     2.387  1.85     9
36   phiBetaT[1, 2]     0.965   1.199    -0.924     1.047     2.654  6.51    15
37   phiBetaT[2, 2]     3.470   1.329     1.301     3.448     5.409  7.89     9
38   phiBetaT[1, 3]     2.774   0.433     1.963     2.773     3.590  1.01    14
39   phiBetaT[2, 3]     6.215   0.670     4.800     6.110     7.637  2.09     8
40   phiBetaT[1, 4]    -0.934   0.372    -1.745    -0.913    -0.228  1.20    19
41   phiBetaT[2, 4]     2.923   0.650     1.698     2.955     4.186  1.01     7
42   phiBetaT[1, 5]    -2.058   0.419    -2.680    -2.132    -1.120  1.97    17
43   phiBetaT[2, 5]     1.407   0.385     0.642     1.411     2.156  2.77    14
44   phiBetaT[1, 6]    -2.382   0.360    -3.038    -2.411    -1.648  1.15    18
45   phiBetaT[2, 6]     1.596   0.388     0.907     1.622     2.377  1.76    19
46   phiBetaT[1, 7]    -0.876   0.610    -2.024    -0.760     0.091  4.48     9
47   phiBetaT[2, 7]     1.836   0.768     0.522     1.678     3.225  1.23     4
48   phiBetaT[1, 8]    -2.552   0.553    -3.638    -2.473    -1.539  2.95    12
49   phiBetaT[2, 8]     0.086   0.303    -0.580     0.122     0.596  1.46    18
50   phiBetaT[1, 9]    -0.128   0.631    -1.379     0.018     1.035  4.05     8
51   phiBetaT[2, 9]     0.674   0.342     0.047     0.674     1.316  1.40    16
52  phiBetaT[1, 10]    -0.032   0.684    -1.677     0.050     0.914  1.07     8
53  phiBetaT[2, 10]     0.003   0.580    -1.079     0.026     1.199  2.86    15
54  phiBetaT[1, 11]     0.605   0.835    -0.869     0.379     2.219  1.07     7
55  phiBetaT[2, 11]    -0.083   0.487    -1.025    -0.059     0.744  1.45    14
56  sumWeightZ01[1]  3722.507  71.368  3579.958  3721.687  3857.958  1.00   574
57  sumWeightZ01[2]  5326.117 124.973  5078.170  5325.332  5576.670  1.23   175
58  sumWeightZ01[3]  7529.181 243.484  7052.642  7538.822  7998.112  1.69    85
59  sumWeightZ01[4] 15466.546 248.563 15010.948 15459.961 15965.160  1.02   219
60  sumWeightZ01[5] 13752.801 328.220 13177.095 13728.943 14476.292  1.01    99
61  sumWeightZ01[6]  9466.853 292.860  8952.612  9437.316 10154.673  1.12    74
62  sumWeightZ01[7]  7608.635 275.206  7101.807  7602.230  8175.634  1.27    69
63  sumWeightZ01[8] 11425.776 528.132 10414.082 11400.303 12494.458  1.32    49
64  sumWeightZ01[9]  4833.239 194.047  4512.097  4812.668  5256.350  1.04   143
65 sumWeightZ01[10]  3107.631 257.927  2718.208  3062.963  3702.399  1.02    42
66 sumWeightZ01[11]  1934.429 186.588  1654.364  1905.694  2395.768  1.60    71
67        sumZ01[1]   983.000   0.000   983.000   983.000   983.000   NaN     0
68        sumZ01[2]  1240.504  18.493  1206.000  1240.000  1279.000  1.44    77
69        sumZ01[3]  1341.041  29.641  1282.975  1341.000  1398.000  1.90    52
70        sumZ01[4]  1408.335  17.593  1376.000  1408.000  1444.025  1.02   185
71        sumZ01[5]  1120.643  20.470  1085.000  1119.000  1165.025  1.01    78
72        sumZ01[6]   672.044  17.758   640.000   671.000   711.000  1.11    52
73        sumZ01[7]   474.720  14.095   447.975   475.000   503.000  1.36    63
74        sumZ01[8]   399.080  16.672   367.000   398.000   433.000  1.37    51
75        sumZ01[9]   221.998   5.219   213.000   221.000   234.000  1.05   131
76       sumZ01[10]    90.636   7.134    80.000    89.000   107.000  1.02    42
77       sumZ01[11]    50.266   4.376    44.000    50.000    61.000  1.53    68

sumWeightZ01ByCohortS <- MCMCsummary(object = d$mcmc, params = c(
        "sumWeightZ01ByCohort"
      ), round = 3) %>%
  rownames_to_column(var = "var")|>
  mutate(
    cohort0 = var |>
      str_extract("\\[(\\d+),") %>%           # Changed to match number before comma
      str_remove_all("\\[|,") %>%             # Remove [ and comma
      as.numeric(),
    ais = var |>
      str_extract(",\\s*(\\d+)\\]") %>%
      str_remove_all(",|\\s|\\]") %>%
      as.numeric(),
    cohort = cohort0 + cohorts[1]
  )

sumZ01ByCohortS <- MCMCsummary(object = d$mcmc, params = c(
        "sumZ01ByCohort"
      ), round = 3) %>%
  rownames_to_column(var = "var") |>
  mutate(
    cohort0 = var |>
      str_extract("\\[(\\d+),") %>%           # Changed to match number before comma
      str_remove_all("\\[|,") %>%             # Remove [ and comma
      as.numeric(),
    ais = var |>
      str_extract(",\\s*(\\d+)\\]") %>%
      str_remove_all(",|\\s|\\]") %>%
      as.numeric(),
    cohort = cohort0 + cohorts[1]
  )

ggplot(sumZ01ByCohortS, aes(ais, mean, group = cohort, color = cohort)) +
  geom_line(aes(linetype = factor(cohort))) +
  labs(#title = "sumZ01 by cohort and age in samples",
       x = "Age in samples",
       y = "Estimated number of fish alive")

ggplot(sumWeightZ01ByCohortS, aes(ais, mean, group = cohort, color = cohort)) +
  geom_line(aes(linetype = factor(cohort))) +
  labs(#title = "sumZ01 by cohort and age in samples",
       x = "Age in samples",
       y = "Estimated sum of fish mass")

19.2 By Cohort, modelNum 2: phi(i,t) * g(i,t), p(i,t) model

Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

      logit(phi[i,t]) <- 
        phiInt[i,t] + 
        phiBeta[1,i,t] * weight[i,t] 

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t)

Survival/growth rate interaction:

Multiplicative

Model runs:

19.2.1 Retrieve model results

library(targets)
library(nimble)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
# 

# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
source('./models/production/modelProduction_OB_functionsToSource.R')
source('./models/production/qmdProduction_OB_functionsToSource.R')

modelNum <- 2 # phi * growth

19.2.2 Model code

# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_2_2002_mostRecent.RData')
d$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:1) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        for (t in 1:T) {
            weightZ01[i, t] <- weightIfAvailable(weight[i, t], 
                sdWeight, meanWeight, z[i, t], available[i, t])
            z01[i, t] <- aliveIfAvailable(z[i, t], available[i, 
                t])
        }
    }
    for (t in 1:T) {
        sumWeightZ01[t] <- sum(weightZ01[1:N, t])
        sumZ01[t] <- sum(z01[1:N, t])
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

getSummaryByCohort <- function(cohort, modelNum = 3) {
  # Load the data
  load(paste0('./models/production/runsOut/production_OB_model', '_', modelNum, '_', cohort, '_mostRecent.RData')) 
  
  # Input data
  eh <- tar_read_raw(paste0('eh_OB_', cohort, '_target'))
  
  print(paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")))

  MCMCplot(object = d$mcmc, params = c("sumZ01"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("sumWeightZ01"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("p"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("grIntT"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("phiIntTOut"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("phiBetaT"), ref_ovl = TRUE)

  get_MCMCTrace_OB(d)
  
   # Create the summary list
  summary_list <- list(
    eh = eh,
    stats = kable(as_tibble(d$runData), caption = "Run statistics"),
    runTime = paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")),
    #traces = get_MCMCTrace_OB(d),
    summaryMod3_tt_growth = MCMCsummary(object = d$mcmc, params = c(
        #"phiIntT", 
        "phiIntTOut", 
        "p", 
        "grIntT", 
        "phiBetaT",
        "sumWeightZ01",
        "sumZ01"
      ), round = 3) %>%
      rownames_to_column(var = "var"),
    summary_OB = get_summary_OB(d, eh)
  )
  
  return(summary_list)
}

#s=summary_OB[, order(summary_OB$ind, summary_OB$t)]

19.2.3 Summaries by cohort

2002

s2002_2 = getSummaryByCohort(2002, modelNum)

[1] "Run time = 21.474 mins"

2003

s2003_2 = getSummaryByCohort(2003, modelNum)

[1] "Run time = 27.403 mins"

2004

s2004_2 = getSummaryByCohort(2004, modelNum)

[1] "Run time = 17.544 mins"

2005

s2005_2 = getSummaryByCohort(2005, modelNum)

[1] "Run time = 8.699 mins"

2006

s2006_2 = getSummaryByCohort(2006, modelNum)

[1] "Run time = 14.192 mins"

2007

s2007_2 = getSummaryByCohort(2007, modelNum)

[1] "Run time = 11.018 mins"

2008

s2008_2 = getSummaryByCohort(2008, modelNum)

[1] "Run time = 9.901 mins"

2009

s2009_2 = getSummaryByCohort(2009, modelNum)

[1] "Run time = 55.498 mins"

2010

s2010_2 = getSummaryByCohort(2010, modelNum)

[1] "Run time = 7.013 mins"

2011

s2011_2 = getSummaryByCohort(2011, modelNum)

[1] "Run time = 6.584 mins"

2012

s2012_2 = getSummaryByCohort(2012, modelNum)

[1] "Run time = 1.09 hours"

2013

s2013_2 = getSummaryByCohort(2013, modelNum)

[1] "Run time = 8.649 mins"

2014

s2014_2 = getSummaryByCohort(2014, modelNum)

[1] "Run time = 6.429 mins"

19.2.4 Combine cohort summaries

sAll_2 <- tibble(.rows = 0) |>
  add_column(!!!s2002_2$summary_OB[0,])

# Fill tibble
for (cohort in cohorts) {
  summary_obj_2 <- get(paste0("s", cohort, "_2"))$summary_OB
  summary_obj_2$cohort <- cohort
  sAll_2 <- bind_rows(sAll_2, summary_obj_2)
}

sAll_2$cohort_ind <- paste0(sAll_2$cohort, "_", sAll_2$ind)

#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all_2 = sAll_2)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

19.2.4.1 Select cohort(s):

summary_OB_OJS_all_T_2 = transpose(summary_OB_OJS_all_2)

import { range } from "d3";

viewof selected_cohort_2 = Inputs.select(d3.range(2002, 2015), {
  label: "Which cohort?",
  value: 2002,
  step: 1,
  multiple: true
})

summary_OB_OJS_all_T_selected_cohort_2 = summary_OB_OJS_all_T_2.filter((d) =>
  (selected_cohort_2.includes(d.cohort))
)

19.2.4.2 Select one or more individuals:

uniqueInds_2 = [...new Set(summary_OB_OJS_all_T_selected_cohort_2.map(d => d.cohort_ind))].sort();

viewof selectInd_mod_2 = Inputs.select(uniqueInds_2, {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})


summary_OB_OJS_all_T_selected_cohort_ind_2 = summary_OB_OJS_all_T_selected_cohort_2.filter((d) =>
  (selectInd_mod_2.includes(d.cohort_ind))
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

19.2.4.3 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Probability of survival" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "pSurv"
      }),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "enc8",
        fill: "orange"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "pSurv"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "first",
        y: 1,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "last",
        y: 1,
        "stroke": "red"
      })
      ,
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2, 
                Plot.selectFirst({
                  x: 11,
                  y: 1,
                  text: d => d.cohort
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_2,
      x: "ind"
    }
  })

Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Standardized body mass (mg)" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "mean_weight"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "mean_weight"
      }),
      <!-- Plot.dot(d, { -->
          <!--   x: "t", -->
          <!--   y: "mean_gr", -->
          <!--   fill: "green" -->
          <!-- }), -->
        Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
          x: "t",
          y: "weightDATA_std",
          fill: "orange"
        }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "first",
        y: 3,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "last",
        y: 3,
        "stroke": "red"
      }),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
                Plot.selectFirst({
                  x: 1,
                  y: 4,
                  frameAnchor: "top-left",
                  text: (d) => "Cohort = " + d.cohort
                })
      ),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
                Plot.selectFirst({
                  x: 1,
                  y: 3.5,
                  frameAnchor: "top-left",
                  text: (d) => "Residual = " + d.meanResid
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_2,
      x: "ind"
    }
  })

19.2.5 Compare cohorts

Combine scohort$summaryMod2_tt_growth data into one data frame

sParamsAll0_2 <- tibble(.rows = 0) |>
  add_column(!!!s2002_2$summaryMod3_tt_growth[0,]) # need to change to mod2

# Fill tibble
for (cohort in cohorts) {
  summary_obj_2 <- get(paste0("s", cohort, "_2"))$summaryMod3_tt_growth # need to change to mod2
  summary_obj_2$cohort <- cohort
  sParamsAll0_2 <- bind_rows(sParamsAll0_2, summary_obj_2)
}

# remove brackets for filtering
sParamsAll0_2$varNoIndex <- sParamsAll0_2$var |>
  str_remove("\\[.*\\]")

# for numeric ordering
sParamsAll1_2 <- sParamsAll0_2 |> filter(varNoIndex != "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+)\\]") %>%
      str_remove_all("\\[|\\]") %>%
      as.numeric()
  )

sParamsAll2_2 <- sParamsAll0_2 |> filter(varNoIndex == "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+),") %>%           # Changed to match number before comma
      str_remove_all("\\[|,") %>%             # Remove [ and comma
      as.numeric(),
    varNumeric2 = var |>
      str_extract(",\\s*(\\d+)\\]") %>%
      str_remove_all(",|\\s|\\]") %>%
      as.numeric()
  )

sParamsAll_2 = bind_rows(sParamsAll1_2, sParamsAll2_2)

19.2.5.1 Combined cohorts

19.2.5.2 p

ggplot(sParamsAll_2 |> filter(varNoIndex == "p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.2.5.3 phi

ggplot(sParamsAll_2 |> filter(varNoIndex == "phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.2.5.4 grIntT

ggplot(sParamsAll_2 |> filter(varNoIndex == "grIntT"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.2.5.5 grIntT, filtered

ggplot(sParamsAll_2 |> filter(varNoIndex == "grIntT", mean > -5, mean < 5), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.2.5.6 phiBetaT

ggplot(sParamsAll_2 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5), aes(varNumeric2, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples") +
  facet_wrap((~varNumeric1))

19.2.5.7 sumWeightZ01

ggplot(sParamsAll_2 |> filter(varNoIndex == "sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.2.5.8 sumZ01

ggplot(sParamsAll_2 |> filter(varNoIndex == "sumZ01"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.2.6 Size-dependent survival

Need to get rid of squared term for cohort-specific models

phiInt_2 <- sParamsAll_2 |> filter(varNoIndex == "phiIntTOut") |>
  dplyr::select(mean, varNumeric1, cohort) |>
  rename(phiIntOut = mean, ais = varNumeric1)

phiBetas_2 <- sParamsAll_2 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5) |>
  dplyr::select(mean, varNumeric1, varNumeric2, cohort) |>
  pivot_wider(names_from = varNumeric1, values_from = mean, names_prefix = "beta") |>
  rename(ais = varNumeric2)

phiSize_2 = expand.grid(weight = seq(-2,2,0.25), ais = 1:11, cohort = cohorts) |> 
   left_join(phiInt_2) |> 
   left_join(phiBetas_2) |>
   mutate(pSurv = ilogit(phiIntOut + beta1*weight))# + beta_2*weight^2))

ggplot(phiSize_2, aes(weight, pSurv, color = factor(cohort), group = cohort)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ais)

19.3 By Cohort, modelNum 3: phi(i,t) * g(i,t)^2, p(i,t) model

Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

      logit(phi[i,t]) <- 
        phiInt[i,t] + 
        phiBeta[1,i,t] * weight[i,t] + 
        phiBeta[2,i,t] * weight[i,t] * weight[i,t] 

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t)

Survival/growth rate interaction:

Multiplicative

Model runs:

19.3.1 Retrieve model results

library(targets)
library(nimble)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
# 

# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
#source('./models/production/modelProduction_OB_functionsToSource.R')
#source('./models/production/qmdProduction_OB_functionsToSource.R')

modelNum <- 3 # phi * growth

19.3.2 Model code

# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_3_2002_mostRecent.RData')
d$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t] + phiBeta[2, i, t] * weight[i, 
                t] * weight[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:2) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
        phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        for (t in 1:T) {
            weightZ01[i, t] <- weightIfAvailable(weight[i, t], 
                sdWeight, meanWeight, z[i, t], available[i, t])
            z01[i, t] <- aliveIfAvailable(z[i, t], available[i, 
                t])
        }
    }
    for (t in 1:T) {
        sumWeightZ01[t] <- sum(weightZ01[1:N, t])
        sumZ01[t] <- sum(z01[1:N, t])
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

19.3.3 Summaries by cohort

2002

s2002_3 = getSummaryByCohort(2002, modelNum)

[1] "Run time = 32.092 mins"

2003

s2003_3 = getSummaryByCohort(2003, modelNum)

[1] "Run time = 37.935 mins"

2004

s2004_3 = getSummaryByCohort(2004, modelNum)

[1] "Run time = 28.262 mins"

2005

s2005_3 = getSummaryByCohort(2005, modelNum)

[1] "Run time = 12.692 mins"

2006

s2006_3 = getSummaryByCohort(2006, modelNum)

[1] "Run time = 20.456 mins"

2007

s2007_3 = getSummaryByCohort(2007, modelNum)

[1] "Run time = 15.911 mins"

2008

s2008_3 = getSummaryByCohort(2008, modelNum)

[1] "Run time = 14.706 mins"

2009

s2009_3 = getSummaryByCohort(2009, modelNum)

[1] "Run time = 1.453 hours"

2010

s2010_3 = getSummaryByCohort(2010, modelNum)

[1] "Run time = 10.1 mins"

2011

s2011_3 = getSummaryByCohort(2011, modelNum)

[1] "Run time = 9.941 mins"

2012

s2012_3 = getSummaryByCohort(2012, modelNum)

[1] "Run time = 1.39 hours"

2013

s2013_3 = getSummaryByCohort(2013, modelNum)

[1] "Run time = 13.827 mins"

2014

s2014_3 = getSummaryByCohort(2014, modelNum)

[1] "Run time = 8.986 mins"

19.3.4 Combine cohort summaries

sAll_3 <- tibble(.rows = 0) |>
  add_column(!!!s2002_3$summary_OB[0,])

# Fill tibble
for (cohort in cohorts) {
  summary_obj_3 <- get(paste0("s", cohort, "_3"))$summary_OB
  summary_obj_3$cohort <- cohort
  sAll_3 <- bind_rows(sAll_3, summary_obj_3)
}

sAll_3$cohort_ind <- paste0(sAll_3$cohort, "_", sAll_3$ind)

#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all_3 = sAll_3)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

19.3.4.1 Select cohort(s):

summary_OB_OJS_all_T_3 = transpose(summary_OB_OJS_all_3)

import { range } from "d3";

viewof selected_cohort_3 = Inputs.select(d3.range(2002, 2015), {
  label: "Which cohort?",
  value: 2002,
  step: 1,
  multiple: true
})

summary_OB_OJS_all_T_selected_cohort_3 = summary_OB_OJS_all_T_3.filter((d) =>
  (selected_cohort_3.includes(d.cohort))
)

19.3.4.2 Select one or more individuals:

uniqueInds_3 = [...new Set(summary_OB_OJS_all_T_selected_cohort_3.map(d => d.cohort_ind))].sort();

viewof selectInd_mod_3 = Inputs.select(uniqueInds_3, {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})


summary_OB_OJS_all_T_selected_cohort_ind_3 = summary_OB_OJS_all_T_selected_cohort_3.filter((d) =>
  (selectInd_mod_3.includes(d.cohort_ind))
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

19.3.4.3 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Probability of survival" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "pSurv"
      }),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "enc8",
        fill: "orange"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "pSurv"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "first",
        y: 1,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "last",
        y: 1,
        "stroke": "red"
      })
      ,
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3, 
                Plot.selectFirst({
                  x: 11,
                  y: 1,
                  text: d => d.cohort
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_3,
      x: "ind"
    }
  })

Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Standardized body mass (mg)" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "mean_weight"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "mean_weight"
      }),
      <!-- Plot.dot(d, { -->
          <!--   x: "t", -->
          <!--   y: "mean_gr", -->
          <!--   fill: "green" -->
          <!-- }), -->
        Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
          x: "t",
          y: "weightDATA_std",
          fill: "orange"
        }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "first",
        y: 3,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "last",
        y: 3,
        "stroke": "red"
      }),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
                Plot.selectFirst({
                  x: 1,
                  y: 4,
                  frameAnchor: "top-left",
                  text: (d) => "Cohort = " + d.cohort
                })
      ),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
                Plot.selectFirst({
                  x: 1,
                  y: 3.5,
                  frameAnchor: "top-left",
                  text: (d) => "Residual = " + d.meanResid
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_3,
      x: "ind"
    }
  })

19.3.5 Compare cohorts

Combine scohort$summaryMod2_tt_growth data into one data frame

sParamsAll0_3 <- tibble(.rows = 0) |>
  add_column(!!!s2002_3$summaryMod3_tt_growth[0,]) 

# Fill tibble
for (cohort in cohorts) {
  summary_obj_3 <- get(paste0("s", cohort, "_3"))$summaryMod3_tt_growth 
  summary_obj_3$cohort <- cohort
  sParamsAll0_3 <- bind_rows(sParamsAll0_3, summary_obj_3)
}

# remove brackets for filtering
sParamsAll0_3$varNoIndex <- sParamsAll0_3$var |>
  str_remove("\\[.*\\]")

# for numeric ordering
sParamsAll1_3 <- sParamsAll0_3 |> filter(varNoIndex != "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+)\\]") %>%
      str_remove_all("\\[|\\]") %>%
      as.numeric()
  )

sParamsAll2_3 <- sParamsAll0_3 |> filter(varNoIndex == "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+),") %>%           # Changed to match number before comma
      str_remove_all("\\[|,") %>%             # Remove [ and comma
      as.numeric(),
    varNumeric2 = var |>
      str_extract(",\\s*(\\d+)\\]") %>%
      str_remove_all(",|\\s|\\]") %>%
      as.numeric()
  )

sParamsAll_3 = bind_rows(sParamsAll1_3, sParamsAll2_3)

19.3.5.1 Combined cohorts

19.3.5.2 p

ggplot(sParamsAll_3 |> filter(varNoIndex == "p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.3.5.3 phi

ggplot(sParamsAll_3 |> filter(varNoIndex == "phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.3.5.4 grIntT

ggplot(sParamsAll_3 |> filter(varNoIndex == "grIntT"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.3.5.5 grIntT, filtered

ggplot(sParamsAll_3 |> filter(varNoIndex == "grIntT", mean > -5, mean < 5), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.3.5.6 phiBetaT

ggplot(sParamsAll_3 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5), aes(varNumeric2, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples") +
  facet_wrap((~varNumeric1))

19.3.5.7 sumWeightZ01

ggplot(sParamsAll_3 |> filter(varNoIndex == "sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.3.5.8 sumZ01

ggplot(sParamsAll_3 |> filter(varNoIndex == "sumZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

19.3.6 Size-dependent survival

Need to get rid of squared term for cohort-specific models

phiInt_3 <- sParamsAll_3 |> filter(varNoIndex == "phiIntTOut") |>
  dplyr::select(mean, varNumeric1, cohort) |>
  rename(phiIntOut = mean, ais = varNumeric1)

phiBetas_3 <- sParamsAll_3 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5) |>
  dplyr::select(mean, varNumeric1, varNumeric2, cohort) |>
  pivot_wider(names_from = varNumeric1, values_from = mean, names_prefix = "beta") |>
  rename(ais = varNumeric2)

phiSize_3 = expand.grid(weight = seq(-2,2,0.25), ais = 1:11, cohort = cohorts) |> 
   left_join(phiInt_3) |> 
   left_join(phiBetas_3) |>
   mutate(pSurv = ilogit(phiIntOut + beta1*weight + beta2*weight^2))

ggplot(phiSize_3, aes(weight, pSurv, color = factor(cohort), group = cohort)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ais)